

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
clearvars -except thetaStar K mpow N reb X Y YY Z specif


%% estimated mean utilities
UStar = payoffs(thetaStar,X,Z,YY,reb,mpow);


%% counterfactuals
load('counterfactualUS_replication.mat','srng_cfactUS')
rng(srng_cfactUS)
NS = 2e+3;
YY4 = YY(YY(:,6)==0,1:5); % 4-mpow cw-intervention strategy profiles
timeUS=0;

% US never intervenes
VnoUS_R = UStar(YY(:,2)==0,1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VnoUS_US = UStar(YY(:,2)==0,2,1:N) + randn(size(YY4,1),NS,N);
VnoUS_UK = UStar(YY(:,2)==0,3,1:N) + randn(size(YY4,1),NS,N);
VnoUS_France = UStar(YY(:,2)==0,4,1:N) + randn(size(YY4,1),NS,N);
VnoUS_Russia = UStar(YY(:,2)==0,5,1:N) + randn(size(YY4,1),NS,N);
VnoUS_China = UStar(YY(:,2)==0,6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
noUSintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameUS',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "UK" "France" "Russia" "China" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VnoUS_UK(2:end,b,n),VnoUS_France(2:end,b,n),VnoUS_Russia(2:end,b,n),VnoUS_China(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_UK=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_France=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_Russia=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_China=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_UK.^YY4(2:end,2)).*((1-sigma_UK).^(1-YY4(2:end,2))).*...
                (sigma_France.^YY4(2:end,3)).*((1-sigma_France).^(1-YY4(2:end,3))).*...
                (sigma_Russia.^YY4(2:end,4)).*((1-sigma_Russia).^(1-YY4(2:end,4))).*...
                (sigma_China.^YY4(2:end,5)).*((1-sigma_China).^(1-YY4(2:end,5)));
            EV_R(m)=VnoUS_R(2:end,b,n)'*EQ(:,m);
        end
        noUSintervEq{b,n}=[EV_R;EQ];
    end
    toc
	timeUS=timeUS+toc;
end
[noUSDMES,noUSintervEq,noUSeqPayoffs] = createDMES(noUSintervEq,VnoUS_R,VnoUS_US,VnoUS_UK,VnoUS_France,VnoUS_Russia,VnoUS_China,[YY4(:,1),zeros(size(YY4,1),1),YY4(:,2:5)],specif);

% US always intervenes
VallUS_R = UStar([1;find(YY(:,2)==1)],1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VallUS_US = UStar([1;find(YY(:,2)==1)],2,1:N) + randn(size(YY4,1),NS,N);
VallUS_UK = UStar([1;find(YY(:,2)==1)],3,1:N) + randn(size(YY4,1),NS,N);
VallUS_France = UStar([1;find(YY(:,2)==1)],4,1:N) + randn(size(YY4,1),NS,N);
VallUS_Russia = UStar([1;find(YY(:,2)==1)],5,1:N) + randn(size(YY4,1),NS,N);
VallUS_China = UStar([1;find(YY(:,2)==1)],6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
allUSintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameUS',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "UK" "France" "Russia" "China" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VallUS_UK(2:end,b,n),VallUS_France(2:end,b,n),VallUS_Russia(2:end,b,n),VallUS_China(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_UK=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_France=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_Russia=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_China=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_UK.^YY4(2:end,2)).*((1-sigma_UK).^(1-YY4(2:end,2))).*...
                (sigma_France.^YY4(2:end,3)).*((1-sigma_France).^(1-YY4(2:end,3))).*...
                (sigma_Russia.^YY4(2:end,4)).*((1-sigma_Russia).^(1-YY4(2:end,4))).*...
                (sigma_China.^YY4(2:end,5)).*((1-sigma_China).^(1-YY4(2:end,5)));
            EV_R(m)=VallUS_R(2:end,b,n)'*EQ(:,m);
        end
        allUSintervEq{b,n}=[EV_R;EQ];
    end
	toc
    timeUS=timeUS+toc;
end
[allUSDMES,allUSintervEq,allUSeqPayoffs] = createDMES(allUSintervEq,VallUS_R,VallUS_US,VallUS_UK,VallUS_France,VallUS_Russia,VallUS_China,[YY4(:,1),ones(size(YY4,1),1),YY4(:,2:5)],specif);


%%
clearvars -except N NS srng_cfactUS YY4 timeUS VnoUS_R VnoUS_US VnoUS_UK VnoUS_France VnoUS_Russia VnoUS_China ...
     VallUS_R VallUS_US VallUS_UK VallUS_France VallUS_Russia VallUS_China noUSintervEq noUSDMES noUSeqPayoffs allUSintervEq allUSDMES allUSeqPayoffs
save('counterfactualUS.mat')



